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The effect of geometric uncertainty due to statistically independent, random, normally 
distributed shape parameters is demonstrated in the computational design of a 3-D 
flexible wing. A first-order second-moment statistical approximation method is used to 
propagate the assumed input uncertainty through coupled Euler CFD aerodynamic / finite 
element stmctural codes for both analysis and sensitivity analysis. First-order sensitivity 
derivatives obtained by automatic differentiation are used in the input uncertainty 
propagation. These propagated uncertainties are then used to perform a robust design of 
a simple 3-D flexible wing at supercritical flow conditions. The effect of the random input 
uncertainties is shown by comparison with conventional deterministic design results. 
Sample results are shown for wing planform, airfoil section, and structural sizing 
variables. 


Introduction 

“A robust design problem is one in which a design is 
sought that is relatively insensitive to uncertain 
quantities” [1]; that is, the objective function does not 
change much and the constraints are not likely to be 
violated. Robust design problems are prevalent in 
disciplines such as structures [2] and mechanisms [3], 
but, due to lack of study to quantify uncertainty in 
aerodynamics and computational expense, have not 
been widely applied for aerodynamic design. This paper 
represents an attempt to study robust design of a 
flexible wing that is built upon an established 
multidisciplinary design optimization procedure. 

Several recent papers [4-8] have addressed the issues of 
and obtained sample results using CFD for 
aerodynamic optimizations subject to uncertainties in 
the input or design variables; we shall call these robust 
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optimizations in contrast to conventional or 
deterministic optimizations. In all of these robust 
optimization demonstrations, three steps are involved. 
First, the input uncertainties are quantified. Second, 
the input uncertainties are propagated through the 
CFD code to obtain uncertainties for the output 
functions. Third, the output functions with 
uncertainties are used in the optimization objective 
and constraint functions to perform a robust design. 
Implementation algorithms and details used to 
accomplish these three steps differ in the cited papers, 
but in all cases the input uncertainty is shown to 
influence the design. These input uncertainties 
represent error sources that are "external" to the CFD 
code. Several other recent papers [9-14] have 
discussed issues and results related to uncertainty 
analysis for CFD codes and applications but have not 
included optimization or design. 

Computational simulation uncertainties also arise from 
physical, mathematical, and numerical modeling 
approximations (see, for example, [11-17]); these are 
called "internal" model error and uncertainty sources. 
Internal sources are not considered herein but must be 
included in assessing the total uncertainty in any 
computational simulation of physical phenomena. The 
recent Journal of Fluids Engineering special section 
on “Quantifying Uncertainty in CFD” [15] contains 
eight invited papers addressing both numerical 
accuracy and physical uncertainty issues. In the 
present work, the discrete code results are considered 
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to be deterministically certain in order to show the 
influence of statistically random input geometry, such 
as may be due to manufacturing variations, on a 
gradient-based performance design optimization. As 
pointed out in [5], [7], and [8], these probabilistic -based 
performance designs differ from the probabilistic-based 
reliability designs that are usually considered in the 
structures disciplines. For the former problem, one is 
more interested in the probable or frequent events that 
occur around the mean values, whereas for the latter 
problem, one is most interested in the infrequently 
occurring catastrophic events. An example of 
multidisciplinary performance optimization subject to 
uncertainty involving linear aerodynamics can be found 
in [18]; simultaneous consideration of both types of 
probabilistic-based design problems is also discussed 
there. 

The present work considers the simultaneous 
aerodynamic-structural performance optimization of a 
simple 3-D flexible wing subject to statistical 
uncertainty in geometric input for wing planform, 
airfoil section, and structural thickness design variables. 
The approach presented in [5] and [8] is extended to 
multidimensional supercritical flow about a flexible 
wing that is modeled by high fidelity multidisciplinary 
computational simulations. A first-order second- 
moment (FOSM) statistical approximation method is 
used to propagate the assumed input uncertainty 
through coupled Euler CFD aerodynamic/finite element 
structural codes for both system analysis and sensitivity 
analysis. First-order sensitivity derivatives (SD) 
obtained by automatic differentiation are used in the 
FOSM input uncertainty propagation. These 
propagated uncertainties are then used to perform a 
robust design of a simple 3-D flexible wing. The effect 
of the random input uncertainties is shown by 
comparing it with conventional deterministic design 
results. This present implementation of the statistical 
approach is very easy to retrofit into gradient-based 
design codes that already utilize analytical or semi- 
analytical sensitivity derivatives for optimization. 

Integrated Statistical Approach 

The integrated statistical approach, as demonstrated for 
a CFD code in [5] and [8] for quasi 1-D subsonic Euler 
flow, is presently implemented for supercritical flow 
about a 3-D flexible wing. This approach follows the 
integrated strategy of [3] that was demonstrated on a 
linkage mechanism design. That strategy for mitigating 
the effect of uncertainty includes (a) uncertainty 
quantification, (b) uncertainty propagation, and (c) 
robust design. Flerein, the details differ somewhat from 
[3], [5], and [8]; the major difference is in the process 
of obtaining the required second-order sensitivity 
derivatives. In this present implementation, finite 
differencing of the "robust" objective and active 


constraint functions is controlled by the optimizer 
code, whereas in the quasi 1-D application [5, 8], these 
second derivative pieces were constructed using both 
hand- and automatically differentiated code as detailed 
in [19] and provided to the optimizer code. 

Uncertainty Quantification 

Uncertainty quantification is not an easy task (for 
example, see [6]); though it is an extremely important 
one that must be included for any realistic design. The 
authors have neither the manufacturing data nor the 
engineering experience to address it here. The 
uncertainty quantification issue was simulated here by 
assuming that the uncertainty in geometry was 
characterized by statistically independent, random, 
normally distributed shape parameters such as those 
due to manufacturing process variations. Our 
assumption simplifies the resulting algebra and 
equations to be coded, but it serves to illustrate the 
influence of input geometry on a 3-D flexible wing 
design. 

Uncertainty Propagation 

Uncertainty propagation was accomplished herein 
using only the FOSM method, whereas in [5] and [8] 
both first- and second-order approximations were 
demonstrated and compared to Monte Carlo 
simulations. This propagation was done for both CFD 
Euler and finite element method (FEM) structures 
discipline codes using the sensitivity analysis codes 
from our Simultaneous Aerodynamic and Structural 
Design Optimization (SASDO) studies [20], First- 
order SDs obtained by automatic differentiation have 
been used as previously demonstrated for a CFD code 
in [5] and [8] and an FEM code in [21], A brief 
outline of the FOSM follows; for more details see [5] 
and [8], 

The first step in the FOSM analyses is to approximate 
the system output solutions of interest in Taylor series 
form. These approximations are formed to estimate the 
output value for small deviations of the input. Given 
input random variables b={b 1; ... ,b n } with means 
b = {b,,...,b n } and standard deviations 
o b = {a bj ,...,a bn } , and system output function F, the 
first-order Taylor series approximations are 

F(b) = F(b) + f^(b j -b i ) (1) 

One then obtains expected values for the mean (first 
moment) and variance (second moment) of the output 
function, F, which depend on the SD and input 
variances <r b . The mean of the output function F , and 
standard deviation O f , are approximated as 
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F = F (b) 
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where the SD, 3F/9 bj , are evaluated at the 
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values b . These SD are obtained using the automatic 
differentiation tool ADIFOR [22, 23], Note that, 
although it is not shown explicitly, the function F may 
also depend on other parameters for which there is no 
uncertainty or for which that uncertainty is too small to 
be significant. Such variables would have no 
contribution to the standard deviation, cr F . They are not 
explicitly included in any further discussion. 


where k is the number of standard deviations a g that 
the constraint g must be displaced to achieve the 
desired or specified target probability P t . For the 
FOSM approximation, the standard deviations ct f and 
are of the form given in Eqs. (2) involving first- 
order SDs. Therefore, a gradient-based optimization 
will then require second-order SDs to compute the 
objective and constraint gradients. Flerein, the 
probabilistic objective and constraint functions of 
Eqs. (4) are evaluated using function analysis and 
first-order SD values from sensitivity analysis codes. 
These augmented functions are then finite differenced 
to obtain the gradients as required by the optimization 
code. 


Robust Design 

Conventional optimization for an objective function 'F, 
which is a function of some system output F, state 
variables Q, and input variables b, is expressed in 
Eqs. (3). 

min x F='F(F;Q,b) 

subject to (3) 

g(F;Q,b) < 0 

The CFD state and structural equilibrium equation 
residuals, R(Q,b)=0, determine Q given b. The system 
constraints g are represented as inequality' constraints. 
The input variables b are precisely known for 
conventional optimization, and all functions of b are 
therefore deterministic. 

For robust design, the conventional optimization must 
be treated in a probabilistic mamier. Given uncertainty 
in the input variables b, all of the functions in Eqs. (3) 
become uncertain. The design variables are now the 
mean values, b = {b,,...,b n } , with b assumed 
statistically independent and normally distributed with 
standard deviations o b . The CFD state and structural 
equilibrium equation residuals R are deemed to be 
satisfied at the mean values Q and b such that 
R(Q,b) = 0 . The objective function is cast in terms of 

expected values and becomes a function of F and cr F . 
The system constraints are cast into a probabilistic 
statement: the probability that the constraints are 
satisfied is greater than or equal to a desired or 
specified target probability, that is P(g < 0) > P t . This 
probability statement is transformed [3] to a constraint 
involving mean values and standard deviations under 
the assumption that the variables involved are normally 
distributed. The robust optimization can be expressed as 
min V F= V F (F, ct f ; Q, b) 
subject to (4) 

gCF;Q,b) + ka ? <0 


Application to 3-D Flexible Wing 

A simple 3-D flexible wing geometry is chosen to 
illustrate the influence that statistically random 
geometric shape and sizing parameters have on a 
gradient-based performance optimization for 
supercritical flow conditions. This application is 
multidisciplinary; it considers both aerodynamics and 
structures. The state vector Q considered above can be 
decomposed into a separate flow-field vector Q and a 
structural displacement vector u, Q={Q,u}. The set of 
design variables b can also be separated into two sets: 
those describing the structural element sizes b slze and 
those describing the wing geometry b ge om such that 
b= !b geom , b size }. The state equations can be separated 
as R-{Ri,j,k,i= R m J where the solution of the high 
fidelitv Euler CFD model is represented by 

Ri.j.k.1 (Q> X def (bgeon. > u))# 0 (5) 

with i, j, k, and 1 representing loop indices over the 
volume mesh points and the mass, momentum, and 
energy equations and the solution of the linear FEM 
structural model is represented by 

R m = K(b size , X jig (b geom ))u - f (Q) = 0 (6) 

with m representing the loop index over the FEM 
mesh points, Xj ig . These equations are coupled in that 
the aerodynamic load f(Q) influences the displacement 
u of the structure and this displacement of the 
aerodynamic shape X de f changes the load. Sensitivity 
analysis essentially involves solving for the derivative 
of Eqs. (5) and (6) with respect to each of the design 
variables b, for dQ/db and5u/Sb . Then the 
derivative equations are coupled by df/db and 
du/db . 

Problem Description 

Two sample problems are considered: a two-design- 
variable (2DV) case for which conventional or 
deterministic optimization results have been 
previously obtained [24], and a four-design-variable 
(4DV) case. The trapezoidal-planform, semispan wing 
and input design variables are shown in Fig. 1. In the 
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2DV examples, the wing sections are held fixed and 
vary linearly from an NACA 0012 at the root to an 
NACA 0008 at the tip, which is then rounded. The size 
of the structural elements is also held fixed in the 2DV 
examples. The two planform design variables are the tip 
setback x t and the tip chord c t . Although these variables 
are not typically subject to significant variance, they are 
ascribed some uncertainty here to illustrate the 
technique and process of propagating input uncertainty 
through the coupled aerodynamic/structural analysis 
and optimizing the design in the presence of the 
uncertainty in the output functions. 



Figure 1. Wing geometry and sizing 
parameterization. 

The relative sizes of the skin thickness, the web 
thicknesses, and the truss cross section areas are fixed 
within each structural zone depicted in Fig. 1. A 
scaling factor F n is assigned to change the thickness and 
area of all structural elements in zone n. In the 4DV 
case, two of these sizing factors are considered as 
design variables. The other two design variables 
considered in the 4DV case were determined by 
examining both sensitivity derivatives and (reasonable) 
expected variances. As can be seen from Eq. (2), the 
size of an output function variance depends on the size 
of the input variances weighted by the sensitivity 
derivative of the output function with respect to that 
input variable. For input variances of 0.1%, the larger 
SDs with respect to root airfoil thickness t r and camber 
z r appeared to provide a more realistic representation of 
practical uncertainty than the variables used in the 2DV 
cases. 

The objective function to be minimized is the negative 
of the square of the lift-to-drag ratio, -(L/D) 2 , for 
M a = 0.8 ,a = 1° , representative of cruise conditions 
for a transport aircraft. Both coupled solution- 


dependent and geometric constraints are imposed. The 
solution-dependent constraints are 

• lower limit on the difference between the total 
lift and the structural weight, L - W , where 
L = C L q„S 

• upper limit on compliance, the work done by the 
aerodynamic loads to deflect the structure, 
V = |J pu ■ n ds 

• upper limit on pitching moment, C m , in lieu of a 
trim constraint 

The purely geometric constraints are 

• minimum leading edge radius, in lieu of a 
manufacturing requirement 

• side constraints (bounds) on the active design 
variables 

State variables for aerodynamics and structures are 
determined from Eqs. (5) and (6), the fluid-flow 
conservation laws (Euler equations) and the structural 
equilibrium conditions (for FEM), respectively. The 
conventional or deterministic optimization for the 
2DV case follows Eqs. (3) and is reported in [24], 
Details regarding procedures, equations, codes, and 
background references are given in that paper and will 
not be repeated here. 

Robust Design Process 

Robust design for both cases follows the integrated 
statistical approach as outlined above, i.e. Eqs. (1), (2), 
and (4), according to [5] and [8], The objective and 
constraint functions listed above are approximated by 
Eqs. (1) with their mean and variances given by 
Eqs. (2). First-order SDs, as obtained from ADIFOR- 
processed analysis code, are used in Eqs. (2) to 
propagate the input uncertainty to the output of this 



Figure 2. Schematic of robust optimization process. 


schematic of this integrated procedure and data flow 
for robust design optimization according to Eqs. (4). 
The shaded portions of the schematic depicted in 
Fig. 2 represent code blocks unchanged from the 
original deterministic gradient-based design code [20, 

24]. 
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In the simple CFD example considered in [5] and [8], 
the objective function was for matching a target, i.e., to 
minimize the square of the difference between an 
output and a predetermined target value. In the 
construction of the objective function obtained by 
taking expected values, the FOSM approximation of the 
variance survives because of the squaring operation. 
Thus, first-order SDs appear in the robust objective 
function. In the present example, the objective function 
is taken as -(L/D) 2 and again the FOSM approximation 
to the variance appears in the robust version. Then 
when the FOSM approximation of the standard 
deviation is used in the constraint, the robust 
optimization problem is given as 


mm -(L/D) 2 


I 

i=l 


d(L/D) ' 


subject to 


( 7 ) 


g = g + k t|Z|^ ff l 


<0 


We note here that the robust optimization problem 
reduces to the conventional or deterministic problem 
when the standard deviation of the inputs CTb is zero. 
Specific computational tools used to perform the tasks 
depicted by the solid boxes in Fig. 2 are identified in 
the next section. 


Computational Tools 

Major computations in this robust optimization 
procedure are performed using a collection of existing 
codes. These codes are executed by a separate driver 
code and scripts. Each code runs independently, some 
simultaneously on separate processors, and the required 
data transfers between them, also directed by the driver, 
are accomplished via data files. 

The aerodynamic flow analysis code used for this study 
is a version of the CFL3D code [25] used in the Euler 
mode. The gradient version of this code, which was 
used for aerodynamic sensitivity analysis, was 
generated by an unconventional application [26] of the 
automatic differentiation code ADIFOR [22, 23] to 
produce a relatively efficient, direct mode, gradient 
analysis code [27], 

The surface geometry was generated based on a code 
utilizing the Rapid Aircraft Parameterization Input 
Design (RAPID) technique developed by Smith et al 

[28] , This code was also preprocessed with ADIFOR to 
generate a code capable of producing SD as well. 

The CFD volume mesh needed by the flow analysis 
code was generated using a version of the CSCMDO 

[29] grid generation code. The associated grid SDs 
needed by the flow sensitivity analysis were generated 
with an automatically differentiated version of 


CSCMDO [30], The 45,000 grid point baseline 
volume mesh required by CSCMDO and used in the 
present flexible wing examples was obtained with the 
Gridgen™ code. The wing surface portion of the mesh 
is shown in Fig. 3. This mesh is admittedly quite 
coarse by current CFD analysis standards. 


CFD mesh 
C-0 topology 
73x25x25 volume 
49x25 on the wing 


FEM mesh 
3251 elements: 

2141 CST 
1110 truss 
583 nodes 

Figure 3. CFD and FEM computational meshes. 

The structural analysis code [3 1] used to compute the 
deflection of the elastic wing was a generic finite 
element code. The flexible structure for the wing 
shown in Fig. 3 was discretized by 583 nodes; there 
were 2,141 constant-strain triangle (CST) elements 
and 1,110 truss elements. Zone boundaries for the 
design variables controlling element size are also 
shown in Fig. 3. Because the elastic deformation was 
assumed to be small, linear elasticity was deemed to 
be appropriate. The structural sensitivity equations 
were derived based on the direct differentiation 
method. The sensitivity of the aerodynamic forces 
appears as a term on the right-hand side (RHS) of the 
deflection sensitivity equations. The derivative of the 
stiffness matrix in these sensitivity equations was also 
generated [32] by using the ADIFOR [22, 23] 
technique. The coefficient matrix of the structural 
sensitivity equations was identical to that of the 
structural equations. Consequently, these structural 
sensitivity equations were solved efficiently by 
backward substitution with different RHSs for each 
sensitivity. 

At the wing surface, i.e., the interface where 
aerodynamic load and structural deflection 

information is transferred, surface nodes of the FEM 
structural model were assumed to be a subset of the 
CFD aerodynamic surface mesh points for this 
application (see Fig. 3). This lack of generality 
allowed for simplifications in the data transfers and, 
although an important issue, it was not deemed crucial 
for these initial 3-D robust optimization 

demonstrations. 

The code for the box labeled Robust System Objective 
and Constraints was simply the programming of Eqs. 
(1), (2) and (4) as previously discussed. These 
probabilistic objective and constraints were then fed to 
the optimizer [33], The Sequential Quadratic 

Programming (SQP) procedure was used and the 
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objective and constraint gradients were calculated by 
the optimizer using finite differences. 

Sample Results & Discussion 

Two-Desian-Variable Cases 

The 2DV optimization problem was solved using 
several values for the input uncertainties associated 
with the design variables. Those results are compared to 
conventional or deterministic design results in Figs. 4. 
Figure 4(a) shows the level sets (contours) of the 
deterministic objective function shaded according to the 
key and drawn as thin solid black lines. In addition, the 
deterministic payload constraint function is shown as a 
yellow (light) solid shaded region to indicate where it is 
violated and as contour lines to indicate where it is 
satisfied. The square symbol is the deterministic 
optimization result. For the deterministic design, the 
payload constraint is active, but the compliance and 
pitching moment constraints are not; in fact, with the 
initial design point of (1,1) they were never active 
during the deterministic optimization process. This 
deterministic result was the initial design point from 
which all the robust optimization problems were 
started. For all the robust design problems, it was 
assumed that the same value ct could be used for each 
of the input c,. 



C, 


Figure 4(a). 2DV level sets of deterministic objective 
and active constraint functions with conventional 
optimization result. 

Figure 4(b) shows robust optimization results overlaid 
on the deterministic results. The robust designs with 
a = 0.01 and a = 0.02 have been pushed away from the 
payload constraint boundary during the optimization 
process. The level sets for the deterministic objective 
function and payload constraint are nearly parallel over 
a significant region so that there is a large allowable 
variation in the design variables to solve the design 
problem to the specified tolerance. This can also be 
recognized by the solution for a second deterministic 
optimization problem shown as the circle for which the 


starting point was the solution of the problem with 
a = 0.02. 



C, 


Figure 4(b). Probabilistic ct = 0.01 & 0.02 robust 
optimization results overlaid on deterministic level 
sets. 

The solutions to the robust optimizations with 
a = 0.03 and a = 0.04 are substantially different from 
the other optimization results as shown in Fig 4(c), 
where they are also overlaid on the deterministic 
results. For these problems, the formerly inactive 
constraint on compliance became active. A red (dark) 
solid region indicates where it is violated and wide red 
(dark) curves indicate where it is not violated. For two 
design variables and two active constraints, the 
solution must be at the intersection of the (robust) 
constraints, which moves the design points much 
farther from the deterministic solution. 



Figure 4(c). Probabilistic a = 0.03 & 0.04 robust 
optimization results overlaid on deterministic level 
sets. 

The effect on the planfonn shape can readily be seen 
in Fig. 5. Although the cr = 0.01 and cr = 0.02 cases are 
easily distinguishable from the deterministic case, they 
are indistinguishable from each other. The a = 0.03 
and ct =0.04 cases are easily distinguishable from all 
other cases. Physically, the effect of reducing the 
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setback is to decrease the twist deformation of the tip 
(wash-out), thereby increasing the lift and hence the 
available payload. But the increased lift, particularly at 
the tip, increases the bending and, with it, the 
compliance. The effect of increasing the tip chord is to 
increase the weight, but that effect is offset by the 
increased lift, which increases the available payload. 
This 2DV problem is not a representative engineering 
problem, but it is useful for understanding the robust 
design process through visualization and also for 
determining that the process is functioning properly. 
Moreover, the size of the coefficients of variation o, 
used here is much larger than one would expect in 
practice for these selected variables. 



Figure 5. Planforms resulting from 2DV 
optimizations. 


Four-Desiqn-Variable Cases 

Although the 2DV problem is useful for understanding 
the process and ascertaining that it is performing 
correctly, the variance of the design variables was 
necessarily exaggerated to see any effect. That is, the 
first step of the robust design process, quantification of 
uncertainty of the input variables, was not followed. For 

Table 1. 4DV problem results. 


this 4DV example, a sensitivity analysis was 
performed to determine those parameters for which a 
realistic variance would have a significant effect on 
the output functions used as the objective and 
constraints. The four parameters chosen as design 
variables were the root airfoil section maximum 
thickness t r , the root airfoil section maximum camber 
Zr, and the structural sizing factors for the two inboard 
regions, and r 2 , as shown in Fig. 1. For the cases 
shown here, a coefficient of variation, CT; = 0.001 , was 
chosen for all four variables. As in the 2DV cases, a 
deterministic optimization was performed first. That 
deterministic result was used as the initial design for 
all the robust optimizations. The deterministic 
optimization process reduced the section thickness t r to 
reduce the shock strength thereby reducing the drag 
and improving the L/D. As a consequence, the wing 
became more flexible. To satisfy the compliance 
constraint g(V), the structure element thickness 
increased, and, as a consequence, the wing became 
heavier. To satisfy the payload constraint g(L-W), the 
section camber increased. The pitching moment did 
not appear to be affected sufficiently for that 
constraint, gfCm), to be active. The deterministic 
results are shown in the first column of Table 1. 

Table 1 and Figs. 6-8 present the results for the 
deterministic optimization and robust optimizations 
for several values of k. An increase of the parameter k 
represents an increase in the specified target 
probability and therefore the probability that the 
constraints are met. Assuming a nonnal Gaussian 
distribution of the output variables, values of k = 1, 2 
and 3 would represent probabilities of 84.13%, 
97.73% and 99.87%, respectively. The mean values of 
the constraint functions are shown in comparison to 
the robust constraint values in Table 1. The design 
variables are compared in Fig. 6. 



Deterministic 

solution 

Robust solution, 
k=1 

Robust solution, 
k=2 

Robust solution, 
k=3 

Desian variable inDuts 




tr 

0.778 

0.774 

0.774 

0.773 

Z r 

1.138 

1.139 

1.139 

1.14 

Ti 

4.075 

4.05 

4.009 

3.984 

r 2 

3.7 

3.656 

3.563 

3.528 

iResDonses 

robust mean 

robust mean 

robust mean 

obj 

-19.1 

-19 

-18.95 

-18.93 

g(L-W) 

-0.00077 

0.000048 -0.0107 

0.00031 -0.0211 

-0.000745 -0.0328 

g(V) 

-0.00032 

-0.000028 -0.0077 

-0.00501 -0.0203 

0.000449 -0.0225 

g(c m ) 

-0.0666 

-0.0531 -0.062 

-0.0422 -0.0601 

-0.0262 -0.0532 

L/D 

21.85 

21.82 

21.77 

21.75 


54243 

53917 

53056 

52684 
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